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Abstract 

The effect of fixed discrete colloidal charges in the primitive model is investi- 
gated for spherical macroions. Instead of considering a central bare charge, as 
it is traditionally done, we distribute discrete charges randomly on the sphere. 
We use molecular dynamics simulations to study this effect on various proper- 
ties such as overcharging, counterion distribution and diffusion. In the vicinity 
of the colloid surface the electrostatic potential may considerably differ from 
the one obtained with a central charge. In the strong Coulomb coupling, 
we showed that the colloidal charge discretization qualitatively influences the 
counterion distribution and leads to a strong colloidal charge-counterion pair 
association. However, we found that charge inversion still persists even if 
strong pair association is observed. 
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I. INTRODUCTION 



The electrostatic interactions in charged colloidal systems play a crucial role in determin- 
ing the physical properties of such materials [|I],0] . The behavior of these systems is extremely 
complex due to the long range Coulomb interactions. A first simplifying assumption is to 
treat the solvent as a dielectric medium solely characterized by its relative permittivity e^. 
A second widely used approximation consists in modeling the short range ion-ion excluded 
volume interaction by hard spheres. These two approximations are the basis of the so-called 
primitive model of electrolyte solutions. The system under consideration is an asymmetrical 
polyelectrolyte made up of highly charged macroions and small counterions in solution. A 
further simplification can be achieved by partitioning the system into subvolumes (cells), 
each containing one macro ion together with its neutralizing counterions plus, if present, ad- 
ditional salt. This approximation has been called accordingly the cell model 0J^. The cells 
assume the symmetry of the macroion, here spherical, and are electrostatically decoupled. 
In this way one has reduced a complicated many-body problem to an effective one colloid 
problem. For spherical macroions the structural charge is normally modeled by a central 
charge, which, by Gauss theorem, is equivalent of considering a uniform surface charge 
density as far as the electric field outside the sphere is concerned. 

Most analytical work as well as simulation approaches rely on the above assumptions. 
It is well known that in the strong Coulomb coupling regime ion-ion correlations become 
very important, and significant deviations from mean-field approaches are expected. One of 
the effects which the mean-field theory like Poisson-Boltzmann can not explain is the phe- 
nomenon of overcharge, also called charge inversion. It consists of binding excess counterions 
to a charged particle (macroion) so that its net charge changes sign. This has recently at- 
tracted significant attention |]5|-[14|. It may give rise to a possible mechanism for strong long 
range attraction between like-sign charged colloids fT^JTSH . 

The purpose of this paper is to investigate if such a phenomenon (overcharge) depends 
on the way the structural charge is represented. The macroion is taken to be perfectly 
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spherical, i. e. we neglect any surface roughness [T^. We introduce discrete charges on 
the macroion sphere instead of a central charge, and compare the results to those obtained 
with a central charge. We concentrate of the following properties in the strong Coulomb 
coupling: overcharging, counterion distribution and surface diffusion. 

II. SIMULATION MODEL 

A. Macroion charge discretization 

The macroion charge discretization is produced by using identical microions of diam- 
eter 0", all identical to the counterions, distributed randomly on the surface of the macroion. 
Then the structural charge is Q = —Z^e = —ZcN^e, where > 0, Z^ is the counte- 
rion valency and e is the positive elementary charge. The discrete colloidal charges (DCC) 
are fixed on the surface of the spherical macroion. In spherical coordinates the elementary 
surface is given by: 

dA = rgsinOdOdLf = —r1d{cos9)d'^ , (1) 

and to produce a random discrete charge distribution on the surface we generated randomly 
the variables cosO and (p. Only configurations leading to an overlap of microions are re- 
jected. Figure |I| shows a schematic view of the setup. Note that in real physical systems 
like sulfonated latex spheres, no large heterogeneities are expected in the charge distribu- 
tion, provided that the colloid surface is relatively regular, therefore our choice is justified. 
Nevertheless, the experimental situation is more complicated since other phenomena such as 
surface chemical reactions fl^, hydration, roughness and many more may be present. 



Here, we restrict ourselves to a simple model in order to understand the effect of macroion 
charge discretization, and leave the other questions for future investigations. 
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B. Molecular dynamics procedure 



We use molecular dynamics (MD) simulations to compute the motion of the counterions 
coupled to a heat bath acting through a weak stochastic force W(t). The equation of motion 
of counterion i reads 

m^ = -V,t/-mr^+W,(t), (2) 

where m is the counterion mass, U is the potential force having two contributions: the 
Coulomb interaction and the excluded volume interaction and F is the friction coefficient. 
Friction and stochastic force are linked by the dissipation-fluctuation theorem < Wi(i) • 
Wj(f) >= 6mrkBTSijS{t — t'). For the ground state simulations the fluctuation force is set 
to zero. 

Excluded volume interactions are taken into account with a pure repulsive Lennard- Jones 
potential given by 



ULj{r) 



4e 



\r-roJ \r-ro J 



+ e, for r - ro < Vcut, , , 

(3) 

0, for r - ro > r^ut-, 



where ro = for the microion-microion interaction (the microion can be a counterion or 
a DCC), ro — 7cr for the macroion-counterion interaction, and Vcut {— is the cutoff 

radius. This leads to a macroion-counterion distance of closest approach a — 8(7. Energy 
and length units in our simulations are defined as e —ksTQ (with To = 298 K), and a — 3.57 
A respectively. 

The pair electrostatic interaction between any pair ij, where i and j denote either a DCC 
or a counterion, reads 

UcUr) = ksToh^ , (4) 
r 

where Ib — / AneoerkBTo is the Bjerrum length describing the electrostatic strength. Being 
essentially interested in the strong Coulomb coupling regime we choose the relative permit- 
tivity tr = 16 {Ib = lOcr), divalent counterions {Zc = 2) and divalent DCC for the remaining 
of this paper. 



The macroion and the counterions are confined in a spherical impenetrable cell of radius 
R. The macroion is held fixed and is located at the center of the cell. The colloid volume 
fraction fm is defined as rf^/R^, where r™, = a — cr/2 is the colloid radius. We have fixed 
R = 40(7 so that fm = 6.6 x 10"'^. To avoid image charge complications, the permittivity er 
is supposed to be identical within the whole cell (including the macroion) as well as outside 
the cell. 



The first step to understand the effect of colloidal charge discretization consists of es- 
timating the electric field, or equivalently, the electrostatic potential generated by such a 
sphere in the absence of counterions. A simple graphical visualization of the field lines is 
here not possible, since there is no perfect symmetry. Indeed, in the present case the electric 
field becomes very anisotropic and irregular close to the sphere, which is the most interest- 
ing region where correlations are expected to be large. To describe qualitatively the effect 
of charge discretization on the electrostatic potential, we compute for three perpendicular 
directions x, y, z the resulting radial potentials Kr(r), Vy{r\ Vz{r) (see Fig. |lD for one given 
DCC random distribution as a function of the distance r > a from the macroion center. 
The radial component of the electric field Ei{r) = — ^V^(r) has the important feature of 
representing the attractive component towards the sphere. The normalized radial potential 
Vi in the i!^^ direction at a distance r from the colloid center is given by: 



where rj{r) is the vector pointing from the microion j to the point where the electric potential 
is computed (see Fig. |1|). Physically, V{r) is the electrostatic potential interaction between 
a counterion and all the surface microscopic colloid charges. The monopole contribution is 



typical bare charges, each corresponding to one given random macroion charge distribution. 



III. MACROION ELECTRIC FIELD 




(5) 



merely given by Vmono{r) 



= —kBTolB ^"!^'' - In Fig. |^ we show the electric potential for three 
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For all cases, one notes that at the vicinity of the surface the potential becomes very different 
from the one computed with a central charge. We carefully checked that similar results 
were obtained for other choices of x, y, z directions (by rotating the trihedron (ex, Cy, Cz)). 
However if we observe the electric field sufficiently far away from the colloidal surface (about 
one macroion diameter) the field is almost exactly the same as the one produced by a central 
charge, which we term isotropic for the rest of this paper. A closer look on Fig. § reveals 
that by increasing the bare charge Zm the electric field starts to become isotropic at smaller 
distances from the sphere's surface. This last feature can be physically easily interpreted. In 
fact when one increases the bare charge, one also increases the absolute number of discrete 
charges which has the effect of approaching the uniform continuous charge density limit 
(corresponding to the isotropic case). 

To capture the discretization effect on the surface electrostatic potential, we have mea- 
sured the electrostatic potential along a circle of radius a concentric to the spherical macroion 
(see Fig. |I]). We start from the top of a given DCC microion and move along a circle in a 
random direction and measure the electrostatic potential V{s) as a function of the arc length 
separation s from the starting point. The same formula as Eq. (^ has been used here. The 
constant monopole contribution is merely given by Vmono = — ^s^o^s"^^^- Results are re- 
ported in Fig. ^for the same configurations as before. It clearly shows that the electrostatic 
potential is strongly fluctuating. More specifically, the higher the structural charge Z^, the 
larger the "oscillation frequency" of the potential fluctuations over the surface. This feature 
can be explained in terms of "/lo/es" . In the very vicinity of a given DCC the potential is 
increased (in absolute value) in average, and around a given DCC there is a hole (depletion 
of charges) which tends to decrease the potential (in absolute value). The average surface of 
this hole is increasing with decreasing bare charge Z^ (i. e. decreasing density of charged 
sites). 

In the following sections we are going to study the effect of charge discretization on the 
counterions distribution in the strong Coulomb coupling. For all following results we used 
the same random charge distributions which gave the results of Figs. ^ and |^. 
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IV. GROUND STATE ANALYSIS 



In this section, we focus on counterion distribution exclusively governed by energy min- 
imization, i. e. T = OK. In such a case correlations are maximal, and all the counterions 
lie on the surface of the spherical macroion. To avoid the trapping in metastable states, 
we systematically heat and cool (10 cycles) the system and retain the lowest energy state 
obtained in this way. Furthermore we choose as the starting configuration one where each 
DCC is exactly associated with one counterion, and each of these dipoles are radially ori- 
ented (each dipole vector and the macroion center lie on the same line). Preliminary, we 
checked that this method reproduces well the ground state energy and structure in simple 
situations where a central charge with two, three, four or five counterions is present. The 
structure of these systems is well known by the Gillespie rules fl^. It turns out that in these 
situations no rough energy landscape (even for = 180 and 90 counterions) appears and 
therefore the MD simulation easily finds the global minimum. It is only in the case of DCC 
that several energy minima are observed. 



First we consider the simple salt free case where the system [macroion + counterions] is 
neutral. In order to characterize the counterion layer, we compute the counterion correlation 
function (denoted by CCF) g{r) on the surface of the sphere, defined as: 



where ps = Nc/'iira'^ is the surface counterion concentration {Nc = Zm/Zc being the number 
of counterions), r corresponds to the arc length on the sphere. Note that at zero temper- 
ature all equilibrium configurations are identical, thus only one is required to obtain the 
CCF. Similarly, one can also define a surface macroion correlation function (MCF) for the 
microions on the surface of the macroion. The CCF is normalized as follows 



A. Neutral case 




(6) 
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/ 2Txrg{r)dr = {N, - 1) . (7) 

JO 

Because of the finite size and the topology of the sphere, g{r) has a cut-off at vra (=25.1 cr). 
Therefore at "large" distance the correlation function differs from the one obtained with an 
infinite planar object. 

The CCF and MCF for two different structural charges (50 and 180) can be inspected 
in Fig. ^. The CCF is computed for a system with a central charge (CC) and for the discrete 
colloid charges (DCC) case. One remarks that both CCF differ considerably following the 
nature of the colloidal charge, i. e., discrete or central (see Fig. For the isotropic 

case (central charge) a Wigner Crystal structure is observed as was already found in Refs. 
||l^ , p!3|JT8| . It turns out that when we have to deal with DCC the counterion distribution is 
strongly dictated by colloidal charge distribution (see Fig. ^). Ground state structures are 
depicted in Fig. ^. It clearly shows the ionic pairing, between DCC and counterions. Also, 
it appears natural to call such a structure a pinned configuration. However one can expect 
that the structure might become less pinned if the typical intra-dipole distance (here a) and 
the typical mean inter-dipole distance become of the same order. This is a case which is not 
discussed in the present paper. It would correspond to extremely highly charged colloids 
that are rarely encountered in nature. Nevertheless, we checked that even for = 360 the 
structure is still pinned, where the average inter-dipole distance is about 2a. 

B. Overcharge 

We now investigate the overcharge phenomenon. The starting configurations corresponds 
to neutral ground states as were previously obtained. The spirit of this study is very similar 
to the one undertaken in Ref. [|1^]. To produce overcharge, one adds successively overcharg- 
ing counterions (OC) at the vicinity of the colloidal surface. Thus the resulting system is 
no longer neutral. By using Wigner crystal theory [|6|,p!9|, we showed that the gain in elec- 
trostatic energy (compared to the neutral state) by overcharging a single uniformly charged 
colloid can be written in the following way [|T^ , |r^ ,p!^ : 
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where AE'^°^', which is equal to the first term of the right member, denotes the gain in energy 
due to ionic correlations for n OC. The functional form of this term can be derived from 
WC theory [p!2| , p!3| , p!8| . The second term on the right hand side, Ai?™°", is the monopole 



repulsion, which sets in when the system is overcharged (with n > 1). This term will, for 
sufficient high number n of OC, stops the process of overcharging. As before Nc = Zm/Zc is 
the number of counterions in the neutral state, and 7 is a constant which was determined 
by using the measured value of AEf'-^ of our simulations. 

The total electrostatic energy of the system as a function of the number of OC is displayed 
Fig. I for four bare charges Z^ (50, 90, 180 and 360 ). These energy curves corresponding 
to discrete systems were produced by averaging over five random DCC realizations. Again, 
the overcharging process is affected by the charge discretization and pinning, but it is still 
energetically favorable. The main effects of charge discretization are: (i) the reduction of 
gain of energy and (ii) the reduction of maximal (critical) acceptance of OC. Both points can 
be explained in terms of ion-dipole interaction. It is exactly this attractive ion-dipole corre- 
lations which is responsible of charge inversion for colloidal systems with discrete charges. 
When the first OC is present, it is normally located in between the pinning centers, and 
will essentially interact with its neighboring dipoles (DCC-counterion). This interaction in- 
creases with decreasing OC-dipole separation, i. e. increasing colloid bare charge Zm- This 
explains why the energy gained increases with Z^ (see Fig. On the other, the repulsion 
between the counterions is not fully minimized since they do not adopt the ideal Wigner 
crystal structure that is obtained with a central charge which in turn explains (i). For a 
higher degree of overcharge, one has to take into account a repulsive monopolar contribution 
identical to AE"^°^ appearing in Eq. (^. Again, since for DCC structures counterions are 
not perfectly ordered, the attractive correlational energy is smaller (in absolute value) than 
the one obtained with a central charge, which in turn explains (ii). Note that for very high 
bare charge {Zm = 360) the overcharge curve obtained with DCC approaches the one from 



the continuous case as expected. 

Common features of overcharging between isotropic and discrete systems are briefly 
given here. We note that the maximal (critical) acceptance of OC (4, 6 and 8 for a central 
charge and 2,4 and 6 for DCC) increases with the macroionic charge (50, 90 and 180 
respectively). Furthermore, for a given number of OC, the gain in energy is always increasing 
with Zra- Also, for a given macroionic charge, the gain in energy between two successive 
overcharged states is decreasing with the number of OC. Note that at T = 0, the value 
acts only as a prefactor. All these features are captured by Eq. (§). 

V. FINITE TEMPERATURE 

In this part, the system is brought to room temperature Tq. We are interested in deter- 
mining the counterions distribution as well as the counterion motion within the counterion 
layer. The radius R of cell is again fixed to 40cr so that the macroion volume fraction 
has the finite value 6.6 x 10~^. Under these conditions the system is still highly energy 
dominated so that at equilibrium all counterions lie on the surface of the macroion (strong 
condensation). 

A. Counterions distribution 

Like in the ground state analysis, we characterize the counterion distribution via its 
surface correlation function. At non zero temperature, correlation functions are computed 
by averaging Z^i^j ^[r — ri)6{T — rf) over 1000 independent equilibrium configurations which 
are statistically uncorrelated. Results are depicted in Fig. |^. For both bare charges Zm (50 
and 180) considered the counterions distribution are affected by charge discretization. The 
effect of temperature is to smooth the CCF. As expected, for the central charge case, the 
counterion positional order is much weaker at room temperature than in the ground state 
case. 
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B. Surface diffusion 



The aim of this section is to answer the following question: do the counterions only 
oscillate around their equilibrium (ground state) position or do they have also a global 
translational motion over the sphere? 

To answer to this question one introduces the following quantity: 

1 

Ax\t,to) = —— dt[x{t)-x{to)]\ (9) 

t — iQ Jto 

which is referred as the mean square displacement (MSD), where x{to) represents the position 
of a given counterion at time t = Iq and a;(t,to) is its position at later time t. The root 
mean square displacement (RMSD) is defined as Aa;(t,to) = y'Ax^(t, to). Like for the 
surface correlation function, the MSD is measured on the spherical surface (arc length) and 
it has a natural upper limit (vra)^. Results for the discrete case are depicted in Fig. || for 
two macroion bare charges (50 and 180), where each counterion' RMSD is given. In 
both cases, the RMSD is smaller than the typical mean inter-dipole separation, which is 
approximatively (ps)~^^^ ~ 6cr for = 50 and (ps)~^^^ ~ 3(t for Z^ = 180. This suggests 
that the motion of the counterion is purely oscillatory around its DCC pinning center. Fig. 
H also shows that pinning is stronger for Z^ = 50 than for Zm = 180. This is in agreement 
with our previous statement, where we point out that the inter-dipole distance has to be 
comparable (or smaller) to (than) the intra-dipole distance in order to reduce pinning effect. 
Thus the DCC sites do produce a noticeable energy well. One can get convinced on this 
point, by evaluating the electrostatic binding energy of an ionic pair Epin which is 

Ep,„ = -kBTohZ^Ja = -AOksn. (10) 

However, for much higher temperatures a liquid like behavior is to be expected. Also, of 
course, the strength of the pinning can be lowered by different parameters: larger ions, higher 
dielectric constant e^, monovalent ions as it is captured by Eq. (|10|). For the isotropic case, 
we have checked that counterions have a large lateral motion and can move all over the 
sphere. This is obvious since in this situation there are no pinning centers. 
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VI. CONCLUDING REMARKS 



We have carried out MD simulations within the framework of the primitive model to 
elucidate the effect of colloidal charge discretization . The main result of our study is that, 
in the strong Coulomb coupling, the charge inversion is still effective when the macroion 
structural charge is carried by discrete charges distributed randomly over its surface area. 
We have shown that the intrinsic electrostatic potential solely due to the surface colloidal 
microions strongly vary from point to point on the macroion sphere. When counterions 
are present, it leads to a pinned structure where every counterion is associated with one 
colloidal charge site. Furthermore we observed a pure phonon-like behavior (counterion 
vibration with small lateral motion) is found at room temperature. 

Future work will address the problem of valency asymmetry, that is when the colloidal 
charges are represented by monovalent counterions and the counterions are divalent. This 
is a non trivial situation since ionic pairing may be frustrated. Also, the case of the low 
Coulomb coupling regime should be addressed. 
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FIGURES 

FIG. 1. Schematic view of the setup: the discrete colloidal charges (DCC) of diameter a are 
in dark grey. The radial electrostatic field components Ex and Ey are represented. For a detailed 
meaning of the other symbols see text. Note that this a a two-dimensional representation of the 
three-dimensional system. 

FIG. 2. Radial electrostatic potential as a function of macroion center distance r produced by 
the fixed microscopic colloidal charges disposed on the sphere. These potential have been measured 
in three perpendicular directions (x, y, z) directions (see Fig. |l|). The isotropic case corresponds 
to the field obtained with a central charge (monopole). Three structural charges are considered: 
(a) Zm = 50 (b) Zm = 90 and (c) Zm = 180. 

FIG. 3. Surface electrostatic potential as a function of the arc length s along a circle of radius 
a concentric to the the macroion for three different trajectories. The monopole contribution is 
represented by the dashed line. The same configurations as those of Fig. |2|(a-c) have been used. 

FIG. 4. Ground state surface correlation functions for two macroion bare charges (a) Z^ = 50 
and (b) Z^ = 180. The two counterion correlation functions (CCF) are obtained for discrete 
colloidal charges (DCC) and for the central charge (CC). To get the same distance range for CCF 
and the colloidal surface discrete microions correlation function (MCF), the MCF curve x-axis 
(r/u) was rescaled by a factor a/ro (compare setup Fig. |l|). 

FIG. 5. Ground state structures for two values (a) Zm = 50 and (b) Z^ = 180 corresponding 
to the two cases of Fig. ^. The colloidal surface microions are in white, and the counterions in 
blue. Full ionic pairing association occurs. 

FIG. 6. Electrostatic energy (in units of fc^To) for ground state configurations of a single 
charged macroion as a function of the number of overcharging counterions for three different bare 
charges Z^- CC stands for the central charge case. The neutral case was chosen as the potential 
energy origin. Dashed lines are produced by using equation (P). 
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FIG. 7. Surface correlation functions at room temperature. The two CCF are obtained for 
discrete colloidal charges (DCC) and for the central charge (CC). (a) Zj^ = 50 (b) Zj^ = 180. 

FIG. 8. Root mean square displacement (rmsd) for each counterion. (a) Zjn = 50 (b) = 180. 
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